1 R Basics

This lecture is to introduce the elements in R programming language that are relevant to decision model and decision analysis.

1.1 Agenda

  • Install R and RStudio
  • Data type
  • Data input and output
  • For loop
  • If statement
  • R function
  • Plots
  • Useful packages in the class
  • Project management

1.2 Install R and RStudio

  • You could install R here.

  • Be aware of your operating system

  • The original R interface is not user-friendly. We recommend using the RStudio IDE here

    • Just download the free version!

  • Open RStudio

1.3 Data type

  • In this class, we will mainly use vector, matrix, array, list, and data.frame.
  • These types of data are used intensively in constructing Markov model and in the package of decision analyses

1.3.1 Vector

x1 <- c(3, 1, 4)
print(x1)
## [1] 3 1 4
x2 <- c(2, 7, 1)
print(x1 + x2)
## [1] 5 8 5
future_human <- c("earther", "martian", "belter")
print(future_human)
## [1] "earther" "martian" "belter"

  • You could combine different type of vectors

    • This operation coerces the data type into text
x1_and_future_human <- c(x1, future_human)
print(x1_and_future_human)
## [1] "3"       "1"       "4"       "earther" "martian" "belter"

  • We often create the initial status in the Markov model using vector

    • The distribution of healthy, sick, and dead in a population
    • You can name the elements in the vector
v_init <- c(0.9, 0.09, 0.01)
names(v_init) <- c("healthy", "sick", "dead")
print(v_init)
## healthy    sick    dead 
##    0.90    0.09    0.01
print(sum(v_init))
## [1] 1
  • Access an element(s) in a vector
print(v_init[2])
## sick 
## 0.09
print(v_init["sick"])
## sick 
## 0.09
print(v_init[2:3])
## sick dead 
## 0.09 0.01
print(v_init[c("sick", "dead")])
## sick dead 
## 0.09 0.01

1.3.2 Matrix & Array

  • Matrix is defined with two dimension using the numbers of columns and rows
x <- matrix(c(1, 0, 0, 
              0.8, 0.2, 0, 
              0, 0, 1), 
            byrow = T, nrow = 3)
print(x)
##      [,1] [,2] [,3]
## [1,]  1.0  0.0    0
## [2,]  0.8  0.2    0
## [3,]  0.0  0.0    1
  • We often use this as the basis of transition matrix in Markov model or use this to trace the outcomes over time.
  • You could also provide the names to the columns and rows
rownames(x) <- colnames(x) <- c("healthy", "sick", "dead")
print(x)
##         healthy sick dead
## healthy     1.0  0.0    0
## sick        0.8  0.2    0
## dead        0.0  0.0    1
  • We can access the value of a specific cell
print(x[2, 2])
## [1] 0.2
print(x["sick", "sick"])
## [1] 0.2

  • Matrix multiplication is often used in Markov models. For example, we can multiply the initial state with the transition matrix

    • v_init is not a matrix. Thus, we use t() to convert v_init into a \(1 \times 3\) matrix.
    • %*% is the symbol for matrix multiplication in R
print(t(v_init) %*% x)
##      healthy  sick dead
## [1,]   0.972 0.018 0.01
  • If transition probabilities in the transition matrix vary with time, array() is more useful.
  • array() can be seen as a data type that stores multiple matrices all at once.
  • Let’s create two transition matrices and combine them into an array.
tr1 <- x
tr2 <- matrix(c(0.9, 0.1, 0, 
                0.7, 0.2, 0.1, 
                0, 0, 1), 
              byrow = T, nrow = 3, 
              dimnames = list(c("healthy", "sick", "dead"), 
                              c("healthy", "sick", "dead")))
  • Create an array with dimension \(3 \times 3 \times 2\)
tr_array <- array(dim = c(3, 3, 2), 
                  data = cbind(tr1, tr2), 
                  dimnames = list(c("healthy", "sick", "dead"), 
                              c("healthy", "sick", "dead"), 
                              c(1, 2)))
print(tr_array)
## , , 1
## 
##         healthy sick dead
## healthy     1.0  0.0    0
## sick        0.8  0.2    0
## dead        0.0  0.0    1
## 
## , , 2
## 
##         healthy sick dead
## healthy     0.9  0.1  0.0
## sick        0.7  0.2  0.1
## dead        0.0  0.0  1.0
  • We can access each transition matrix by proper indexing
print(tr_array[ , ,  1])
##         healthy sick dead
## healthy     1.0  0.0    0
## sick        0.8  0.2    0
## dead        0.0  0.0    1
print(tr_array[ , ,  2])
##         healthy sick dead
## healthy     0.9  0.1  0.0
## sick        0.7  0.2  0.1
## dead        0.0  0.0  1.0
  • We can also access a specific value of the transition array
print(tr_array[1, 2, 1])
## [1] 0
  • Similarly, we can perform calculation using a slice of array
v_time2 <- t(v_init) %*% tr_array[ , , 1]
print(v_time2)
##      healthy  sick dead
## [1,]   0.972 0.018 0.01
v_time3 <- v_time2 %*% tr_array[ , , 2]
print(v_time3)
##      healthy   sick   dead
## [1,]  0.8874 0.1008 0.0118

1.3.3 List

  • List can combine different types of data without affecting the data type.
temp_ls <- list(future_human = future_human, v_init = v_init, tr_array = tr_array)
print(temp_ls)
## $future_human
## [1] "earther" "martian" "belter" 
## 
## $v_init
## healthy    sick    dead 
##    0.90    0.09    0.01 
## 
## $tr_array
## , , 1
## 
##         healthy sick dead
## healthy     1.0  0.0    0
## sick        0.8  0.2    0
## dead        0.0  0.0    1
## 
## , , 2
## 
##         healthy sick dead
## healthy     0.9  0.1  0.0
## sick        0.7  0.2  0.1
## dead        0.0  0.0  1.0
  • There are two approaches to access an element in the list
print(temp_ls[[1]])
## [1] "earther" "martian" "belter"
print(temp_ls$v_init)
## healthy    sick    dead 
##    0.90    0.09    0.01
  • You can access a specific element in an element of a list
print(temp_ls$tr_array[, , 1])
##         healthy sick dead
## healthy     1.0  0.0    0
## sick        0.8  0.2    0
## dead        0.0  0.0    1

1.3.4 data.frame

  • If your Stats or Biostats classes use R, you probably have encountered data.frame() very frequently.
  • Let’s create a data.frame()
profile <- data.frame(name = c("Amos", "Bobbie", "Naomi"), 
                      human_type = c("earther", "martian", "belter"), 
                      height = c(1.8, 2.1, 1.78))
print(profile)
##     name human_type height
## 1   Amos    earther   1.80
## 2 Bobbie    martian   2.10
## 3  Naomi     belter   1.78
  • data.frame() is essentially a named list of vectors with the same length.
typeof(profile)
## [1] "list"
length(profile$name) 
## [1] 3
length(profile$human_type)
## [1] 3
length(profile$height)
## [1] 3
  • Access elements in data.frame()
profile[profile$name == "Bobbie", ]
##     name human_type height
## 2 Bobbie    martian    2.1
profile[profile$name == "Bobbie", "height"]
## [1] 2.1

1.4 Path to your working directory

  • It is recommended to create a folder/directory for each of your project.
  • You would store all your R script, data, and outputs in the directory.
  • Whenever you are working on a project, set up your working directory first.
setwd("path-to-your-project-folder")
setwd("zoeslaptop/Documents/CEA/introR/")
  • The setwd() function sets up your working directory and allows you to access the files and folders in the working directory much easier.

  • For example, I want to run the R script hello.R in the Rscript/ folder in your working directory.

    • Without setting up the working directory first, I have to type source("zoeslaptop/Documents/CEA/introR/Rscript/hello.R")
    • If setting up the working directory first, I can run the script by this command: source("Rscript/hello.R")

  • The format of path could be slightly different between windows, unix, and linux

  • RStudio provides a very nice feature project. A project sets up the working environment automatically. Whenever you open a project, you are in the working directory to a specific project alreadly!

  • You could find more information about creating an R project and other useful tip here.

1.5 Data input and output

  • You could load/read or save/write your own datasets.

  • R data format family:

  • .RData, .rda, and .rds.
  • You can retrieve the information in these data objects using functions load() and readRDS()

  • You can save the data using save() and saveRDS().

  • It is common that you might encounter other data types (e.g., .csv and .txt) or even data format for other software (e.g., Stata and SAS).

  • To read or save other types of data, you could find the information here and here.

1.6 Loops

  • Loops help repeat routine or calculations.
  • Let’s watch a video from Codecademy.

  • Components included in loops

    • A well-defined routine/process
    • When to start
    • When to stop

  • The most used loop is for loop in many decision models. The structure follows:

for (i in c(start : end)) {
  # Routine / Process
}

Example 1: Fibonacci sequence is the sum of the two preceding numbers. We start the sequence from 0 and 1. What are the first 10 Fibonacci numbers? (0 and 1 are the 1st and 2nd Fibonacci number, respectively.)

  • Note that we want to get all the 10 Fibonacci numbers. We need to create a vector to store all the 10 numbers.

    • We already know the first two values. Thus, we replace the first two values in the fib_vec.
fib_vec <- rep(0, 10)
fib_vec[c(1:2)] <- c(0, 1)
print(fib_vec)
##  [1] 0 1 0 0 0 0 0 0 0 0
  • The start point and end point are 3 and 10 because we already knew the first two Fibonacci numbers.
for (i in c(3 : 10)) {
  fib_vec[i] <- fib_vec[i - 1] + fib_vec[i - 2]
}

print(fib_vec)
##  [1]  0  1  1  2  3  5  8 13 21 34

Example 2: In year 2400, there are 3000 martians in Mars colony. The growth rate of the martian population follows this formula \(0.05 P(t) \bigl(1 - \frac{P(t)}{10000})\). What is the population size in year 2500?

  • Note that the growth rate depends on the population size at each year.
  • You can use for loop to calculate the population size at each year.
# initialize a vector of population size over the next 100 years
popsize <- rep(3000, 100)

# Calculate 
for (t in c(1 : 100)) {
  popsize[t + 1] <- popsize[t] + 0.05 * popsize[t] * (1 - popsize[t] / 100000)
}

popdf <- data.frame(year = c(2400:2500), 
                    popsize = popsize)
print(popdf[popdf$year == 2500, ])
##     year  popsize
## 101 2500 81499.86

1.7 If statement

  • If statement is used when you just want to execute the chunk of the code if a certain condition is met.
  • The structure follows
if (condition1) {
  # Execute some code 
} else if (condition2) {
  # Execute some code
} else {
  # Execute some code
}
  • else if and else are not always needed.
  • Before we jump into examples of if statement, let’s take a look at how to write the condition in the if()
"yoda" == "windu"
## [1] FALSE
Jedi <- c("yoda", "windu", "kenobi")
"yoda" %in% Jedi
## [1] TRUE
  • If statement only execute the code when the condition in if() is true.
Mandalorian <- c("satine", "sabine", "jango")

x <- "yoda"

if (x %in% Jedi) {
  print("May the force be with you!")
} else if (x %in% Mandalorian) {
  print("This is the way!")
} else {
  print("Hello, world!")
}
## [1] "May the force be with you!"
  • Note that R is case sensitive. If you type "Yoda", you get "Hello, world!"

1.8 R function

  • As your code grows, copying the entire decision model many times is annoying.
  • Consider programming some routines / processes in a function if you will reuse the routines / processes many times.
  • There are many example functions in R, e.g., lm(), sum(), print(), etc.

  • The primary components of an R function are the body(), formals(), and environment(). We often want to return the results from the function.

    • body(): the code inside the function.
    • formals(): the input arguments.
    • environment(): where the variables in the function are located.
    • return(): return the relevant output
speak <- function(x) {
  Jedi <- c("yoda", "windu", "kenobi")
  Mandalorian <- c("satine", "sabine", "jango")
  
  membership <- "Not Jedi or Mandalorian"
  if (x %in% Jedi) {
    say <- "May the force be with you!"
    membership <- "Jedi"
  } else if (x %in% Mandalorian) {
    say <- "This is the way!"
    membership <- "Mandalorian"
  } else {
    say <- "Hello, world!"
  }
  return(list(say = say, membership = membership))
}

formals(speak)
## $x
body(speak)
## {
##     Jedi <- c("yoda", "windu", "kenobi")
##     Mandalorian <- c("satine", "sabine", "jango")
##     membership <- "Not Jedi or Mandalorian"
##     if (x %in% Jedi) {
##         say <- "May the force be with you!"
##         membership <- "Jedi"
##     }
##     else if (x %in% Mandalorian) {
##         say <- "This is the way!"
##         membership <- "Mandalorian"
##     }
##     else {
##         say <- "Hello, world!"
##     }
##     return(list(say = say, membership = membership))
## }
environment(speak)
## <environment: R_GlobalEnv>
  • What is the result return from the function?
speak("jango")
## $say
## [1] "This is the way!"
## 
## $membership
## [1] "Mandalorian"
  • In this class, we might ask you to change some part of the code in a markov model function.

Question: How would you program the Matian population growth in a function? What are the input arguments? What are the results return from the function?

1.9 Plots

1.9.0.1 Histogram

  • Showing the distribution of a variable
x <- rnorm(1000, 0, 1) # draw 1000 samples from a normal distribution
print(mean(x))
## [1] -0.02976383
print(sd(x))
## [1] 0.996775
hist(x, freq = F, col = "gray")

1.9.0.2 Scatter plots

  • Inspect the correlation between two variables
library(CEAutil)
data(worldHE)

print(head(worldHE))
##        Entity Code Year LEyr HealthExp     Pop
## 1 Afghanistan  AFG 1800   NA        NA 3280000
## 2 Afghanistan  AFG 1801   NA        NA 3280000
## 3 Afghanistan  AFG 1802   NA        NA 3280000
## 4 Afghanistan  AFG 1803   NA        NA 3280000
## 5 Afghanistan  AFG 1804   NA        NA 3280000
## 6 Afghanistan  AFG 1805   NA        NA 3280000
  • Only focus on the data in 2010
worldHE2010 <- worldHE[worldHE$Year == 2010, ]
hist(worldHE2010$LEyr, main = "Life expectancy in 2010")

hist(worldHE2010$HealthExp, main = "Health expenditure in 2010 per capita")

plot(worldHE2010$HealthExp, worldHE2010$LEyr, ylim = c(70, 90), 
     xlab = "health expenditure", ylab = "life expectancy")
text(worldHE2010$HealthExp, worldHE2010$LEyr, labels = worldHE2010$Entity, cex = 0.8)

source

1.9.0.3 Line plots

  • What is the trend of health expenditure in the US 1970-2015?
worldHE_US <- worldHE[worldHE$Entity == "United States" & worldHE$Year >= 1970 & worldHE$Year <= 2015, ]
# adding lines: UK's expenditure at the same time period
worldHE_UK <- worldHE[worldHE$Entity == "United Kingdom" & worldHE$Year >= 1970 & worldHE$Year <= 2015, ]

plot(worldHE_US$Year, worldHE_US$HealthExp, type = "l", col = "red", 
     ylab = "expenditure", xlab = "year", ylim = c(500, 10000))
lines(worldHE_UK$Year, worldHE_UK$HealthExp, col = "royalblue")
legend("topleft", legend=c("US", "UK"),
       col=c("red", "royalblue"), lty = 1, cex = 0.8)

  • There is a great package to help you make better plots! ggplot2

1.10 Useful packages in the class

  • The packages that we will use in the class include: ggplot2, dampack, and CEAutil.
  • Here are the code to install these packages.
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(devtools)) install.packages("devtools")
if(!require(remotes)) install.packages("remotes")
if(!require(dampack)) remotes::install_github("DARTH-git/dampack", dependencies = TRUE)
if(!require(CEAutil)) remotes::install_github("syzoekao/CEAutil", dependencies = TRUE)

2 Decision tree

2.1 Agenda

  • Introduction of Amua
  • Creating a decision tree using Amua
  • Example
  • Exporting Amua decision tree to R script
  • Prepare the R decision tree from Amua for further analyses

For a simple decision tree, we can draw the tree easily. As the decision tree grows, there is software helping you draw the tree such as TreeAge and Amua. In this class, we use Amua because it is free. We will show how to use Amua to build a decision tree and export the tree to R script for CEA analysis.

2.2 Install Amua

Follow the instruction here to install Amua. Be aware of the difference between Mac and Windows users! After you install Amua, remember where you store the software.

2.3 Create a decision tree

2.3.1 Open the Amua.jar.

2.3.2 Decision, chance, terminal nodes, and decision tree

  • Amua includes three type of nodes: decision, chance and terminal nodes in a decision tree.

2.3.3 Parameterization

  • While you can hard code all the parameter values in each branch or terminal node, you could parameterize the input parameters (especially the parameters that might vary).

2.3.4 Adding more outcomes

  • You could add multiple outcomes of interest at the terminal nodes.
  • Sometimes we might be interested not only in cost but also in epidemiological outcomes

2.3.5 Change analysis

  • Change your analysis type to cost-effectiveness analysis.
  • Remember to select the cost and effectiveness outcomes.
  • Also select the baseline strategy.

2.3.6 Save the model and export the tree to R code

  • We then save the model and export the tree structure to R code.
  • Notice that Amua created a folder, called temp_Export/.
  • In this folder, there are two .R files: main.R and functions.R.
  • The body of the decision tree model is in main.R.

2.4 Example

Dracula is preparing for a spring break party tonight. But there’s one final decision that Dracula is struggling with. Being a vampire, he intends to bite and suck the blood of some of his guests at the party (about 25% of them, by his best estimate). While being bitten by a vampire won’t turn his guests into vampires or zombies, like some horror movies might suggest, there is a 50% chance that a vampire bite results in a rather severe bacterial infection, of which 66% of cases require hospitalization for an average of 1 night.

Being a gracious host, Dracula is considering different ways of administering antibiotic prophylaxis to his guests to reduce their risk of infection. One option is to administer the antibiotics to his victim‐guests just before he bites them – this would reduce the risk of a vampire bite infection by 20%. However, the antibiotics can be even more effective, reducing the risk of infection by 90%, if administered at least 30 minutes before being bitten. To achieve this, Dracula is considering putting the antibiotics into the drinks served at the party to ensure that his guests are all properly dosed before he bites his victims. But this means that all his guests would be exposed to the antibiotics (not just those he intends to bite), and he knows that about 5% of people are severely allergic to these antibiotics and would require immediate hospitalization if exposed. Dracula is therefore also considering not administering antibiotic prophylaxis at all to avoid this harm.

However, all the healthy blood doesn’t come without cost. This party will cost Dracula $1000 for a total of 200 guests (an average of $5 per guest). In addition, Dracula expects the cost of antibiotic to be $10 for each guest he bites. If Dracula decides to administer the antibiotics in all the drinks served at the party, the total cost of antibioitics is expected to be $100 (Dracula gets discount for buying a large batch of antibiotics!). Also, Dracula is willing to pay for the cost of hospitalization for any guest who experience bacteiral infection due to his bites because he feels responsible. The cost of hospitalization per person per night is $500. Dracula expects to get an average of 470 mL healthy blood by biting a guest. However, if the guest ends up hospitalized due to bacterial infection, the healthy blood that Dracula can get is reduced by 10%. Dracula hopes that you can help him determine which strategy he should choose.

Think about the following quesitons:

  1. List out Dracula’s strategies.
  2. What are the potential outcomes for Dracula’s party?
  3. Draw the decision tree for Dracula’s party.

In-class activities:

  1. How to convert decision tree from Amua to R code and calculate outcomes.
  2. How to vary parameter values and obtain outcomes for different parameter values.